#ifndef __MATH_COMPLEX_H__
#define __MATH_COMPLEX_H__

/*
===============================================================================

Complex number

===============================================================================
*/

namespace inf {
    class Complex {
    public:
        float				r;		// real part
        float				i;		// imaginary part

        Complex( void );
        Complex( const float r, const float i );

        void 				set( const float r, const float i );
        void				zero( void );

        float				operator[]( int index ) const;
        float &				operator[]( int index );

        Complex			operator-() const;
        Complex &			operator=( const Complex &a );

        Complex			operator*( const Complex &a ) const;
        Complex			operator/( const Complex &a ) const;
        Complex			operator+( const Complex &a ) const;
        Complex			operator-( const Complex &a ) const;

        Complex &			operator*=( const Complex &a );
        Complex &			operator/=( const Complex &a );
        Complex &			operator+=( const Complex &a );
        Complex &			operator-=( const Complex &a );

        Complex			operator*( const float a ) const;
        Complex			operator/( const float a ) const;
        Complex			operator+( const float a ) const;
        Complex			operator-( const float a ) const;

        Complex &			operator*=( const float a );
        Complex &			operator/=( const float a );
        Complex &			operator+=( const float a );
        Complex &			operator-=( const float a );

        friend Complex	operator*( const float a, const Complex &b );
        friend Complex	operator/( const float a, const Complex &b );
        friend Complex	operator+( const float a, const Complex &b );
        friend Complex	operator-( const float a, const Complex &b );

        bool				compare( const Complex &a ) const;						// exact compare, no epsilon
        bool				compare( const Complex &a, const float epsilon ) const;	// compare with epsilon
        bool				operator==(	const Complex &a ) const;						// exact compare, no epsilon
        bool				operator!=(	const Complex &a ) const;						// exact compare, no epsilon

        Complex			    reciprocal( void ) const;
        Complex			    sqrt( void ) const;
        float				abs( void ) const;

        int					get_dimension( void ) const;

        const float *		to_float_ptr( void ) const;
        float *				to_float_ptr( void );
        const char *		to_string( int precision = 2 ) const;

        static Complex ORIGIN;
    };

    INF_INLINE Complex::Complex( void ) {
    }

    INF_INLINE Complex::Complex( const float r, const float i ) {
        this->r = r;
        this->i = i;
    }

    INF_INLINE void Complex::set( const float r, const float i ) {
        this->r = r;
        this->i = i;
    }

    INF_INLINE void Complex::zero( void ) {
        r = i = 0.0f;
    }

    INF_INLINE float Complex::operator[]( int index ) const {
        assert( index >= 0 && index < 2 );
        return ( &r )[ index ];
    }

    INF_INLINE float& Complex::operator[]( int index ) {
        assert( index >= 0 && index < 2 );
        return ( &r )[ index ];
    }

    INF_INLINE Complex Complex::operator-( void ) const {
        return Complex( -r, -i );
    }

    INF_INLINE Complex &Complex::operator=( const Complex &a ) {
        r = a.r;
        i = a.i;
        return *this;
    }

    INF_INLINE Complex Complex::operator*( const Complex &a ) const {
        return Complex( r * a.r - i * a.i, i * a.r + r * a.i );
    }

    INF_INLINE Complex Complex::operator/( const Complex &a ) const {
        float s, t;
        if ( math::fabs( a.r ) >= math::fabs( a.i ) ) {
            s = a.i / a.r;
            t = 1.0f / ( a.r + s * a.i );
            return Complex( ( r + s * i ) * t, ( i - s * r ) * t );
        } else {
            s = a.r / a.i;
            t = 1.0f / ( s * a.r + a.i );
            return Complex( ( r * s + i ) * t, ( i * s - r ) * t );
        }
    }

    INF_INLINE Complex Complex::operator+( const Complex &a ) const {
        return Complex( r + a.r, i + a.i );
    }

    INF_INLINE Complex Complex::operator-( const Complex &a ) const {
        return Complex( r - a.r, i - a.i );
    }

    INF_INLINE Complex &Complex::operator*=( const Complex &a ) {
        *this = Complex( r * a.r - i * a.i, i * a.r + r * a.i );
        return *this;
    }

    INF_INLINE Complex &Complex::operator/=( const Complex &a ) {
        float s, t;
        if ( math::fabs( a.r ) >= math::fabs( a.i ) ) {
            s = a.i / a.r;
            t = 1.0f / ( a.r + s * a.i );
            *this = Complex( ( r + s * i ) * t, ( i - s * r ) * t );
        } else {
            s = a.r / a.i;
            t = 1.0f / ( s * a.r + a.i );
            *this = Complex( ( r * s + i ) * t, ( i * s - r ) * t );
        }
        return *this;
    }

    INF_INLINE Complex &Complex::operator+=( const Complex &a ) {
        r += a.r;
        i += a.i;
        return *this;
    }

    INF_INLINE Complex &Complex::operator-=( const Complex &a ) {
        r -= a.r;
        i -= a.i;
        return *this;
    }

    INF_INLINE Complex Complex::operator*( const float a ) const {
        return Complex( r * a, i * a );
    }

    INF_INLINE Complex Complex::operator/( const float a ) const {
        float s = 1.0f / a;
        return Complex( r * s, i * s );
    }

    INF_INLINE Complex Complex::operator+( const float a ) const {
        return Complex( r + a, i );
    }

    INF_INLINE Complex Complex::operator-( const float a ) const {
        return Complex( r - a, i );
    }

    INF_INLINE Complex &Complex::operator*=( const float a ) {
        r *= a;
        i *= a;
        return *this;
    }

    INF_INLINE Complex &Complex::operator/=( const float a ) {
        float s = 1.0f / a;
        r *= s;
        i *= s;
        return *this;
    }

    INF_INLINE Complex &Complex::operator+=( const float a ) {
        r += a;
        return *this;
    }

    INF_INLINE Complex &Complex::operator-=( const float a ) {
        r -= a;
        return *this;
    }

    INF_INLINE Complex operator*( const float a, const Complex &b ) {
        return Complex( a * b.r, a * b.i );
    }

    INF_INLINE Complex operator/( const float a, const Complex &b ) {
        float s, t;
        if ( math::fabs( b.r ) >= math::fabs( b.i ) ) {
            s = b.i / b.r;
            t = a / ( b.r + s * b.i );
            return Complex( t, - s * t );
        } else {
            s = b.r / b.i;
            t = a / ( s * b.r + b.i );
            return Complex( s * t, - t );
        }
    }

    INF_INLINE Complex operator+( const float a, const Complex &b ) {
        return Complex( a + b.r, b.i );
    }

    INF_INLINE Complex operator-( const float a, const Complex &b ) {
        return Complex( a - b.r, -b.i );
    }

    INF_INLINE Complex Complex::reciprocal( void ) const {
        float s, t;
        if ( math::fabs( r ) >= math::fabs( i ) ) {
            s = i / r;
            t = 1.0f / ( r + s * i );
            return Complex( t, - s * t );
        } else {
            s = r / i;
            t = 1.0f / ( s * r + i );
            return Complex( s * t, - t );
        }
    }

    INF_INLINE Complex Complex::sqrt( void ) const {
        float x, y, w;

        if ( r == 0.0f && i == 0.0f ) {
            return Complex( 0.0f, 0.0f );
        }
        x = math::fabs( r );
        y = math::fabs( i );
        if ( x >= y ) {
            w = y / x;
            w = math::sqrt( x ) * math::sqrt( 0.5f * ( 1.0f + math::sqrt( 1.0f + w * w ) ) );
        } else {
            w = x / y;
            w = math::sqrt( y ) * math::sqrt( 0.5f * ( w + math::sqrt( 1.0f + w * w ) ) );
        }
        if ( w == 0.0f ) {
            return Complex( 0.0f, 0.0f );
        }
        if ( r >= 0.0f ) {
            return Complex( w, 0.5f * i / w );
        } else {
            return Complex( 0.5f * y / w, ( i >= 0.0f ) ? w : -w );
        }
    }

    INF_INLINE float Complex::abs( void ) const {
        float x, y, t;
        x = math::fabs( r );
        y = math::fabs( i );
        if ( x == 0.0f ) {
            return y;
        } else if ( y == 0.0f ) {
            return x;
        } else if ( x > y ) {
            t = y / x;
            return x * math::sqrt( 1.0f + t * t );
        } else {
            t = x / y;
            return y * math::sqrt( 1.0f + t * t );
        }
    }

    INF_INLINE bool Complex::compare( const Complex &a ) const {
        return ( ( r == a.r ) && ( i == a.i ) );
    }

    INF_INLINE bool Complex::compare( const Complex &a, const float epsilon ) const {
        if ( math::fabs( r - a.r ) > epsilon ) {
            return false;
        }
        if ( math::fabs( i - a.i ) > epsilon ) {
            return false;
        }
        return true;
    }

    INF_INLINE bool Complex::operator==( const Complex &a ) const {
        return compare( a );
    }

    INF_INLINE bool Complex::operator!=( const Complex &a ) const {
        return !compare( a );
    }

    INF_INLINE int Complex::get_dimension( void ) const {
        return 2;
    }

    INF_INLINE const float *Complex::to_float_ptr( void ) const {
        return &r;
    }

    INF_INLINE float *Complex::to_float_ptr( void ) {
        return &r;
    }


}   // ::inf

#endif /* !__MATH_COMPLEX_H__ */
